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We have proposed new algorithms for the numerical integration of the equations of motion 
for classical spin systems. In close analogy to symplectic integrators for Hamiltonian equations of 
motion used in Molecular Dynamics these algorithms are based on the Suzuki- Trotter decomposition 
of exponential operators and unlike more commonly used algorithms exactly conserve spin length 
and, in special cases, energy. Using higher order decompositions we investigate integration schemes 
of up to fourth order and compare them to a well established fourth order predictor-corrector 
method. We demonstrate that these methods can be used with much larger time steps than the 
predictor-corrector method and thus may lead to a substantial speedup of computer simulations of 
the dynamical behavior of magnetic materials. 
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I. INTRODUCTION 

Collective phenomena in many materials can be traced back to the presence of interacting magnetic moments on 
the atomic level. The exploration of magnetic systems therefore plays a key role in both physics and materials science. 
The understanding of phase transitions, critical phenomena, and scaling has in part been founded on the investigation 
of simple model spin systems such as the Ising, the XY, and the Heisenberg model (see below). These spin systems 
continue to be of high relevance for the investigation of dynamic critical behavior and dynamic scaling. On the other 
hand realistic models of magnetic materials can be constructed from these simple spin models, if the interaction 
parameters and the underlying lattice structure are taken as an input from experiments. However, in such cases 
the theoretical analysis of experimentally accessible quantities, such as the dynamic structure factor, is usually too 
demanding for analytical methods. Computer simulations of the dynamical behavior of spin systems have, therefore, 
become a very important tool for the theoretical understanding of dynamic critical behavior and material properties 
of magnetic systems, as the following examples may show. 

Large scale computer simulations have been performed in recent years in order to explore the dynamical behavior of 
classical XY jj| and Heisenberg models [|| in d — 2 dimensions and of classical Heisenberg ferro- and antiferromagnets 
in d — 3 HQ]. The simulations are based on model Hamiltonians for continuous degrees of freedom represented by a 
three-component spin with fixed length |S&| = 1 for each lattice site k. A typical model Hamiltonian is then given 

by 

n = J E ( s k s f + s i s i + xs kSn - D ^2(Sk) 2 , (i.i) 

<k,l> k 

where J is the exchange integral, < k,l > denotes a nearest- neighbor pair of spins Sfc, A is an exchange anisotropy 
parameter, and D determines the strength of a single-site or crystal field anisotropy. For A = 1 and D — Eq.([0]) 
represents the classical isotropic Heisenberg ferromagnet or the corresponding antiferromagnet for J > or J < 0, 



respectively. In the case A = D = Eq.(l.l) reduces to the XY model. 
. Realistic descriptions of specific magnetic materials may require additional interactions in the Hamiltonian, like 
an additional two-spin exchange interaction between next nearest neighbors, third nearest neighbors, etc. (see, e.g., 
Ref. |^]). Two spin interactions do not always provide a sufficient representation of the interactions in a magnetic 
■ system. A more accurate description of the isotropic ferromagnet EuS for example also requires a three spin exchange 
interaction Q of the type 

n 3 = -j 3 Yl (Sfc ■ so (s fe • s m ) , (i.2) 

<fc,/,m> 

where < k, I, m > denotes a triple of nearest neighbor spins. Four spin exchange coupling constants are often negligible, 
but the biquadratic interaction given by 
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H 4 = -J 4 (Sfc'S/) 2 , (1-3) 



<k,l> 



which is well known from the Blumc - Emery - Griffith (BEG) model (j, needs to be considered in certain cases |J. 

The thermodynamic properties of th e mo del can be obtained from a Monte-Carlo simulation of the Hamiltonian 
7~Ltot — T~i + ^3 + T~t4 given by Eqs. (|l.l| ), (1.2), and (fL3|). In order to study the dynamic properties of the spin system 



the equations of motion given by [tlH4 



di Sk = "asT x s * (L4) 



must be integrated nume rica lly, where a Monte-Carlo simulation of the model provides equilibrium configura tion s as 



initial conditions for Eq.( |l.4| ). Note that frequencies will be measured in energy units so that K = 1 in Eq.(L4). A 
typical quantity to be determined from the dynamics of the model is the dynamic structure factor <S(q, cj), which is 
given by the space-time Fourier transform of the spin-spin correlation function 

g^{r k -T h t-t>) = {S%{t)SP{t')), (1.5) 

where a, (3 = x,y, z denote the spin component, and r/ are lattice vectors, and the average (. . .) must be taken over a 
sufficiently large number of independent initial equilibrium configurations. The effect of collective thermal excitations 
of the system (e.g., phonons) is not included in the above approach. However, in magnetic systems the typical time 
scale on which Eq.( |1.4| ) has to be integrated is much shorter than typical time scales set by these excitations, so that 
the above procedure can be justified. In general one has to resort to the solution of Langevin equations rather than 



Eq.(1.4) in order to obtain the correct dynamics, but this is beyond the scope of this article. 

From the above outline of a typical Monte-Carlo spin dynamics study it is evident, that at least in d — 3 the 
maj or p art of the CPU time needed for a spin dynamics simulation is consumed by the numerical integration of 
Eq.( |l.4[ ). It is therefore desirable to do the integration with the biggest possible time step. However, using standard 
methods a severe restriction on the size of the time step is posed by the accu racy within which the numerical method 
observes the conservation laws of the dynamics. It is evident from Eq.( |l.4[ ) that |Sfc| for each lattice site k and the 
total energy are conserved. Symmetries of the Hamiltonian impose additional conservation laws. If, for example, 



'Htot — T~t according to Eq.(l.l) for D = and A = 1 (isotropic Heisenberg model) the magnetization M = ^ fc Sfc is 
a conserved quantity. For the anisotropic Heisenberg model, i.e., A ^ 1 and/or / only the z-component M z of 
the magnetization is conserved. Conservation of spin length and energy is particularly crucial, because the condition 
| Sfc | = 1 is a major part of the definition of the model and the energy of a configuration determines its statistical 
weight. It would therefore also be desirable to devise an algorithm which conserves these two quantities exactly. 

The outline for the remainder of this investigation is as follows. In Sec. II we describe the properties of the currently 
used fourth-order predictor-corrector method Sec. Ill is devoted to the presentation of our new integration 

procedure, based on Trotter-Suzuki decompositions of exponential operators, which conserve spin length and, for the 
isotropic case only, energy exactly. In Sec. IV we compare both schemes with special regard to the accuracy within 
which the conservation laws hold and their speed for a given demand on accuracy. A comparison within the full 
framework of a spin dynamics simulation for a Heisenberg ferromagnet is also presented. We finally discuss the 
prospects of both algorithms in view of large scale spin dynamics simulations in Sec.V. 

II. PREDICTOR-CORRECTOR METHODS 

Predictor-corrector methods provide a very general tool for the numerical integration of initial value problems 
like Eq.( |1.4| ) with an equilibrium configuration as the initial value. The demand for very good spin length and 
energy conservation, however, requires methods with small truncation errors in the time step St. Further stability 
requirements for long times have led to the implementation of a fourth-order scheme, which we briefly reproduce here 
for the convenience of the reader. 



In a more symbolic form Eq.(L4) can be written as y = f(y) with the initial condition y(0) = yo, where y is a 
short-hand notation of a complete spin configuration written, e.g., as a 3iV dimensional vector for a lattice with N 
sites. The initial equilibrium configuration is simply denoted by yo- The predictor step of the scheme is then given, 
e.g., by the explicit Adams-Bashforth four-step method Q 

y(t + St) = y(t) + ^ [55/(y(t)) - 59f(y(t - St)) + 37/(j/(f - 2St)) - 9f(y(t - 3St))} (2.1) 
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which has a local truncation error of the order (St) 5 . 
implicit Adams-Moulton three-step method || 

St 



The corrector step consists of typically one iteration of the 



y(t + St) = y(t) + — [9f(y(t + St)) + 19 f(y(t)) - 5f(y(t - St)) + f(y(t - 2St))} 



(2.2) 



which also has a local truncation error of the order (St) 5 
knowledge of y(St), y (2St) , and y(3St) apart from y(0) = y 
of y 



The first application of Eq.(2.1) obviously requires the 
These can be provided by tl 



iree successive integrations 

f(y) (see Eq.(L4)) by the fourth-order Runge-Kutta method || which cannot be used for the entire time 
integration, because truncation err ors accu mula te too fast during typical integration times required for high frequency 
resolution. In order to apply Eqs.(f2~l|) and (2^2) repeatedly the spin configuration at the last four time steps must be 
kept in memory. 

It is interesting to note that the predictor-corrector method according to Eqs. (|2.l|) and ( |2.2| ) observes the conserva- 
tion laws for the magnetization according to the symmetry of the Hamiltonian to within machine ac cura cy if peri odic 
boundary conditions are employed. The reason is that all spins are updated simultaneously by Eqs.(2.1) and Eq.(2.2) 
and that the driving forces (torques) are evaluated precisely according to the exact equation of motion (see Eq.(1.4)) 
for every time step. The periodic boundary conditions restore the discrete translational invariance of the infinite 
lattice in the finite sample needed for the simulation. 

The predictor-corrector method is very general and, therefore, its implementation is independent of the special 
structure of the right-hand side of the equation of motion. Apart from the conservation of the magnetization the 
other conservation laws discussed in Sec. I will therefore only be observed within the accuracy set by the truncation 
error of the method. In pr actic e, this limits the time step to typically St = 0.01/ J in d = 3 j| for the isotropic model 
Hamiltonian given by Eq.(l.l) (D = 0), where the total integration time is of the order 200/ J. The same method 
has also been used in d = 2 with a time step St = 0.025/ J g, but in this case the total integration time did not 
exceed 60/ J. In view of the fact that a time res olut ion of typically 0.1/ J — 0.2/ J is desired for the evaluation of the 
time displaced correlation function given by Eq.(1.5), it is clear that the total CPU time required for a spin dynamics 
simulation in d = 3 is dominated by the numerical integration of the equations of motion. A quantitative analysis is 
presented in Sec. IV. 



III. SUZUKI-TROTTER DECOMPOSITION METHODS 



The motion of a spin according to Eq.(1.4) may be visualized as a Larmor precession of the spin S around an 
effective axis il which is itself time d epe ndent. In order to illustrate our reasoning we restrict ourselves again to the 
Hamiltonian Htot = T~t giv en b y Eq.(l.l) for the simple case D = but for arbitrary values of A. The evaluation of 
the right-hand side of Eq. ( |l.4| ) then shows that the lattice can be decomposed into two sublattices such that a spin 
on one sublattice performs a Larmor precessi on i n a local field Q of neighbor spins which are all located on the other 
sublattice. For the Hamiltonian given by Eq.(Ll) there are only two such sublattices if the underlying lattice is, e.g., 
simple cubic or bcc. 

The first basic idea of the algorithm to be described below is to actually perform a rotation of a spin about its 
local field f2 by an angle a — \fl\St, rather than integrating Eq.( |1.4| ) by some standard method. This procedure 
guarantees the conservation of the spin length |S| to within machine accuracy. (A procedure like this has already been 
implemented in a numerical study o f sp in motions in a classical ferromagnet [[L0j). The second basic idea is to exploit 
the sublattice decomposition of Eq.(1.4) to also ensure energy conservation to within machine accuracy. Denoting the 
two sublattices by A and £>, respectively, we can write Eq.(1.4) in the form 



dt 



SkeA = ^e[{S}] x S ke A 



dt 



S keB = £!^[{S}] x S fce g, 



(3.1) 



where J7^[{S}] and f2g[j_Sj] denote the local fields produced by the spins on sublattice A and £>, respectively. Either 
of the equations in Eq^J^l) reduces to a linear system of differential equations if the spins on the other sublattice are 
kept fixed. This suggests an alternating update scheme, i.e., the spins S^g^i are rotated for the given values of SkeB 
and vice versa. This implies that the scalar products S^g^ • f2e[{S}] remain constant during the update of Sfcg^ and 
the scalar products S^gg • r2^[{S}] remain constant during the update of S fcee . From Eq.(l.l) for D = and Eq.(1.4) 
we therefore find that the energy is exactly conserved during this alternating update scheme (see also Eq^^)). Note, 
that each sublattice rotation is performed with the actual values of the spins on the other sublattice, so that only a 
single copy of the spin configuration is kept in memory at any time. However, the magnetization will not be conserved 
during the above rotation operations. Moreover, the two alternating rotation operations do not commute, so that a 
closer examination of the sublattice decomposition of the spin rotation is required. 
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In order to obtain a simple description of the operations performed on the spin configuration during the integration 
of the equations of motion we again represent a full configuration by a vector y wh ich is decomposed into two sublattice 
components yj( and ys according to y = (y^, ys)- The cross products in Eq.(3.1) can be expressed by matrices A and 
B which are the infinitesimal generators of the rotation of the spin configuration y A on sublattice A at fixed y& and 
of the spin configuration ys on sublattice B at fixed y^, respectively. The exact update of the configuration y from 
time t to t + St can then be expressed by an exponential (matrix) operator according to 



y(t + St) = e^ st y(t). 



(3.2) 



Although the exponential operator in Eq.(3.2) rotates each spin of the configurati on it has no simple explicit form, 
because the rotation axis for each spin depends on the configuration itself (see Eq.(3.1)) and is therefore not known 



a priori. However, the operators e Adt and e BSt which rotate y_4 at fixed yg and ys at fixed y_4, respectively, do have 



a simple explicit form. We demonstrate this for the case A = 1 and D = in Eq.(l.l). For each k G A we find 



n A [{s}] = -j = 

l=NN(k) 



(3.3) 



where NN(k) denotes the nearest neighbors of k ( whic h belo ng t o yg). Eq^^) can be readily generalized for A ^ 1, 
the case D ^ will be discussed below. From Eqs.(pTj) and (3.3) we obtain (see also Rcf. p0[|) 



S k (t + St) 



flk(^k ■ Sfc(t)) 



s fc (t) 



fi fc (n fc -s fc (t)) 



cos(\Qk\5t) 



Ofc x S fc ft) 



sin(|n fc |Jt). 



(3.4) 



Note, that acc ording to Eq.(3.4) 0,% ■ Sfc(i + St) = ■ Sfc(i) which explicitly confirms energy conservation. Fo r 
k € B EqJ3.4) also holds in exactly the same form. The alternating updat e sch eme for the integration of Eq.( |3.f|) , 
i.e., Eq.(jl.4|) now amounts to the replacement e ( A + B ) St — > e A8t £ BSt nl Eg. ( |3.2| ) , which is only correct up to terms 
of the order (St) 2 pj| . The magnetization will therefore only be conserved up to terms of the order St (global 
truncation error), which is insufficient for practical purposes. The remedy for this shortcoming, which constitutes 
the third basic idea of the algorithm, is now obvious. We simply employ higher order Suzuki- Trotter decompositions 
of the exponential operator in Eq.(3.2) to increase the order of the local truncation error of the algorithm and thus 
improve the magnetization conservation. The simplest possible improvement is given by the well known second order 
decomposition [Tl] ] 



,(A+B)5t 



e ASt/2 e BSt e A8t/2 



0{St 3 



(3.5) 



which will be used for comparison with the predictor-corrector method outlined in Sec. II. Note, that Eq.( |3.5| ) is 
equivalent to the midpoint integration method applied to Eq.(BJ) (see also Ref. [|l(|). We furthermore use the fourth 
order decomposition III] 



e (A+B)St = Y^ e p l ASt/2 e p z BSt e p i A5t/2 + 



(3.6) 



with the parameters 



Pi = P2 = Pi = Pb = P = 1/(4 



and P3 = 1 — 4p. 



(3.7) 



It is also possible to construct a decomposition like Eq . (pi] ) wit h on ly three factors, but then [p,-| > I for i = 
1,2,3 which is numerically unfavorable |Tl|. From Eqs. Q3.5| ) and (3.6) the close analogy to symplectic integrators 
obtained from the Liouville operator formalism for Molecular Dynamics simulations is obvious [ |f 2| . The sublattice 
decomposition of the spin degrees of freedom in spin dynamics corresponds to the natural decomposition of the degrees 
of freedom in Hamiltonian dynamics into positions and momenta. Especially the second order decomposition given 
by Eq.([3.5|) is equivalent to the velocity Verlet algorithm for Molecular Dynamics simulations |T^] (see also Ref. |l^] 
for recent developments in this field) which itself is equivalent to the well-known leapfrog algorithm |T^ |. 



As will be sho wn in Sec. IV , th e addi tion al computational effort to be invested in the evaluation of Eq.(3.6) as 
compared to Eq.(3.5) or Eqs.(2T) and ( |2.2|) can be compensated by using larger time steps. The evaluation of the 
trigonometric functions in Eq.(3.4) can also be avoided by observing that the above decompositions are only correct to 
within a certain order in St. It is therefore sufficient to replace sin a: by its Taylor polynomial T(x), where T(x) = x if 
Eq.(3.5) is used and T(x) — x — x 3 /6 if Eq.(3.6) is used. The cosine in Eq.(3.4) then has to be replaced by yl — T 2 (x) 



4 



in order to maintain conservation of spin length. It is also worth noting that the above decompositions maintain the 
time inversion property of e^ A+BS>St , i.e., their inverse is obtained by the replacement St — > —St. A Hamiltonian with 
both nearest and next-nearest neighbor two-spin interactions on, e.g., simple cubic or bcc lattices, can be treated 
within the above framework if the lattice is decomposed into four sublattices. In contrast to the predictor - corrector 
method conservation of magnetization according to the symmetry of th e Ha milto nian is only observed within the 
truncation error of the decomposition method, because according to Eqs.( |3.5|) and (3.6) the sublattices A and B are 
no longer strictly equivalent. 

We now generalize the above considerations to the case D ^ in Eq.fluj). For a spin in sublattice A the equation 
of motion reads 



dl 



S kf zA — f2g[{S}] x S keA - 2DS keA e z x S keA 



(3.8) 



where denotes the unit vector pointing alon g the z-axis, and spins in sublattice B obey an equation of the same form. 
In contrast to the isotropic case (see Eq.(3.1)) the equation of motion for each individ ual spin on each sublattice is 
nonlinear. It is possible to obtain an analytic solution of this equation, in the spirit of Eq.( [3~4| ), where the trigonometric 
functions are basically replaced by Jacobian elliptic functions, and to implement this solution as a rotation operation 
in one of the above decomposition schemes. This approach, however, involves very precision sensitive operations which 
greatly limit speed and therefore the overall efficiency of such an algorithm. In practice, it is much more convenient 
to include the effects of the nonlinearity in Eq.(3.8) by an iterative method, as described in the following paragraph. 



For the sublattice decomposition of the spin rotation in the isotropic case discussed above the requirement for 
energy conservation in the presence of a single site anisotropy reads 



n k ■ s k (t + st)-D [s z k (t + st)f = n k • s fc (t) - d [s z k {tf 



(3.9) 

for k S A and k £ B, where fl k is given by Eq. fl3.3p . In order to perform a rotation operation in analogy to 
Eq.(3.4) we have to identify an effective rotation axis. This can be achieved by rewriting Eq.( |3.9| ) in the form 
fife • (Sfe(i + St) — Sfe(t)) = 0, where Q k is given by 



fife = 0* - D (0, 0, S z k (t) + S z k (t + St)) , 



(3.10) 



i.e., in order to perform the rotation S k at the future time t + St must be kno wn in advance. This problem can only 
be solved iteratively starting from the initial v alue S k (t + S t) — S k (t) in Eq.( 3.10 ) and performing several updates 
according to the decompositions given by Eqs.(3.5) or fl3.6| ), respectively, in order to improve energy conservation 
according to Eq.(3.9). Naturally, this procedure leads to a substantial slowdown of the integration algorithm, where 
the ene rgy is no longer exactly conserved. Details are discussed in Sec. IV. The biquadrati c in teraction given by 
Eq.(L3) can be treated by the same iterative scheme. The three-spin interaction given by Eq.([l.2|), however, requires 
a reconsideration of the sublattice decomposition. Specific numerical examples in comparison with those obtained 
from the predictor-corrector method are presented in the following section. 



IV. THE METHODS IN COMPARISON 



For a quantitative analysis of the integration methods outlined above we restrict ourselves to the Hamiltonian 
T~ltot = 7i given by Eq.([0]) for A = 1 in d = 3. The underlying lattice is simple cubic with L = 10 lattice sites in 
each direction and periodic boundary conditions are imposed in all cases discussed below. The simulations have been 
performed on IBM RS/6000 workstations at the Center for Simulational Physics and a Cray C90 at the Pittsburgh 
Supercomputer Center. 

In order to compare the different integration methods we first investigate the accuracy within which the conservation 
laws are fulfilled. As a standard initial c onfi guration we choose a well equilibrated configuration from a Monte-Carlo 
simulation of the model defined by Eq.(Ll) for A = 1 at a temperature T — 0.8T C for D = and D = J, where 
T c refers to the critical temperature of the isotropic model (D = 0) and is given by K c = J/(fcgT c ) = 0.693035 
Ifufl . The magnetization of such a configuration differs from zero and provides an indicator for the numerical quality 
of the magnetization conservation. We integrate the equations of motion to t = 800/ J and monitor the energy 
e(t) = E(t)/(JL 3 ) of the configuration per spin and in units of the exchange coupling constant J and the modulus 
m(t) = |M(t)|/L 3 of the magnetization per spin for the isotropic case D = and its z-component m z (t) = M z (t)/L 3 
for the stongly anisotropic case D — J as functions of time. Note, that for these tests both integration methods are 
started from identical initial configurations. The time step for the predictor-corrector method is chosen as St — 0.01/ J 
in all cases. 
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We first consider the case of D = for which the implementation of Eqs. fl3.5|) and ( |3.6| ) using Eq.( |3.4| ) is straight- 
forward. The Taylor polynomial T(x) is chosen as T(x) = x — x 3 /6 for Eq.( |3.5| ) and T(x) = x — x 3 /6 + x 5 /120 for 
Eq.( |3.6[ ) in order to reduce the amplitude of the magnetization fluctuations. Figjl] shows that e(i) for the predictor- 
corrector method increases linearly with time whereas the decomposition methods both yield a constant value for e(t) 
as shown in Fig.|l|. Fig.|] displays the magnetization conservation for the predictor-corrector method (St = 0.01/ J) 
and the second order decomposition method according to Eq. (|3.5| ) for St = 0.04/ J. The predictor-corrector method 
conserves m(t) exactly, whereas the second order decomposition causes fluctuations of m(t) on all time scales. The 
apparent deviation of the temporal average of the magnetization from the exact value is due to the finite sample 
time; for times greater than 2000/ J the deviation changes sign and fluctuations on time s cale s larger than 800/ J 
bec ome visible. The temporal structure of m(t) for the decomposition methods given by Eq.( |3.5| ) for St = 0.04/ J and 
Eq.Q for St = 0.2/ J is displayed in Fig|. It is remarkable that the fourth order decomposition gives a substantially 
better magnetization conservation than the second order decomposition despite the very large time step. In order to 
achieve the same overall accuracy of the magnetization conservation with the secon d or der decomposition a time step 
St < 0.02/ J is neccessary. A single integration of the equations of motion using Eq.([3~5|) (second order decomposition) 
is a bout twice as fast as the predictor-corrector method used here. The more complex fourth order decomposition (see 
Eq.(3.6)), however, is about 2.5 times slower than the predictor corrector method. Taking the increase in time step by 
factors of 4 and 20, respectively, into account, we find that both decomposition methods yield an eightfold speedup 
of the integration of the equation of motion. If the overall quality of the magnetization cons ervation is also taken 
into account, there is a clear advantage for the fourth order decomposition according to Eq.(3.6) for the isotropic case 
D = 0. 

We now turn to the strongly anisotropic case D = J. The predictor-corrector method can be applied as before, 
but the decomposition scheme needs a maj or modification because the spin rotation axis depends on the spin value 
S% at the future time t + St (see Eg. ( 3. 10] )). As already pointed out in Sec. Ill this gives rise to a self consistency 
problem which can be solved iteratively, where the quality of th e en ergy conservation depends on the number of 
iterations performed. For the second order decomposition (see Eq.(|3.5|)) and the time step St — 0.04/ J two iterations 
are sufficient to obtain a better energy conservation than the predictor-corrector method. This means that a single 
integration of the equations of motion using the second order decomposition takes twice as long as for D = 0, so that 
its advantage in speed only comes fro m t he increase in the time step, which is still a factor of four. For the fourth 
order decomposition according to Eq.(3.6) with St = 0.2/ J six iterations are needed to obtain energy conservation to 
within six significant digits, which makes one integration with the fourth order decomposition 15 times slower than 
one integration with the predictor-corrector method. From the increase of St by a factor 20 only a 30% gain in speed 
remains for Eq.( |3.6| ). The number of iterations needed decreases with St, but this decrease does not compensate 
the loss in speed due to the smaller time step. However, one still obtains a greatly improved energy conservation. 
The direct comparison of the energy conservation is shown in Fig.|[ All three methods behave in a similar way, 
the change in energy is basically linear with time. The reason for this is the iterative nature of all three methods 
in the case D ^ 0. The magnetization fluctuations shown in Fig.|| for the predictor-corrector and the second order 
decomposition method (St = 0.04/ J) reveal the same behavior of both methods as shown in FigJH for D = 0. The 
direct comparison between the second and fourth order decomposition for St — 0.04/ J and St — 0.2/ J, respectively, 
is displayed in Fig.[| The overall accuracy of the magnetization conservation appears to be independent of D for 
both decomposition metho ds. If the emphasis is put on overall energy conservation and speed, the second order 
decomposition given by Eq.( |3.5|) has some advantages over the predictor-corrector method. If the emphasis is put on 
energy conservation alone, the fourth order decomposition given by Eq.([T^) shows the best performance, but it is 
only slightly faster than the predictor-corrector method. 

We close this section with a direct comparison of the longitudinal (5;(q, ui)) and tra nsverse (St(q, oj)) components 
of the dynamic structure factor for the isotropic Heisenberg ferromagnet (see Eq.(l.l) for D = 0) on a simple cubic 
lattice (L = 10, periodic boundary conditions) at the temperature T = 0.8T C in the (100) direction. In order to 
measure iS;(q, u>) one has to consider the projections of all spins onto the direction of the magnetization M of the 
initial configuration. Note that in general M for a given initial (equilibrium) spin configuration will be nonzero even 
in a finite system. Only the average (M) of the magnetization over sufficiently many equilibrium configurations yields 
(M) = within the statistical errors. The spin components transverse to M then give access to ${(q, oS). The results 
are shown in Figs.|7j and || for q = (7r/5, 0, 0), where 5;(q, to) and St(q, w) have been normalized to the static structure 
factor <Sz,t(q) = J £;,t(q, u>)du. The diamonds represent the result for the predictor-co rrec tor method (St — 0.01/ J) 
and the triangles show the result for the second order decomposition method (see Eq.( |3.5| ), St = 0.04/ J). For both 
methods the equations of motion have been integrated to 800/ J and averages have been taken over 1000 initial 
configurations, where the time displaced correlation functions have been measured to 400/J. The statistical error 
indicated by the error bar represents one standard deviation. For the longitudinal structure factor shown in Fig]?] 
the statistical error is smaller than the symbol size. The overall agreement of the data is very good. The main 
maximum corresponds to the spin wave peak and is located at cjq = 0.25 J in all cases. The frequency resolution is 
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5u> = 2ir J/400 ~ 0.016 J. The shoulder like feature at w ~ 0.5 J in <S/(q, ui) (see Fig.^) is due to many-spin wave 
processes, the description of which is beyond the scope of this article. Both methods have spent about the same CPU 
time on generating the initial configurations, but with the time steps given above, the integration of the equations 
of motion with the second order decomposition method is about eight times faster than with the predictor-corrector 
method. The discrepancy between the errorbars shown in Fig.|| is of statistical origin and indicates the statistical 
uncertainty of the errorbars themselves. Note, however, that the discrepancy only occurs in a frequency region well 
outside the spin wave peak, where the signal level is almost three orders of magnitude below its peak value. 

V. SUMMARY 

Starting from a well developed and established predictor-corrector method for the numerical integration of the 
equations of motion of a classical spin system as a timing and precision standard we devised and tested an alternative 
class of algorithms which is based on Suzuki- Trotter decompositions of exponential operators. Our findings can be 
summarized as follows. 

1. Advantages of the predictor-corrector method 

The greatest advantage of the predictor-corrector method is its versatility and its capability to conserve the 
magnetization exactly. Isotropic and anisotropic spins systems with nearest and next-nearest neighbor inter- 
actions can be treated within one and the same numerical approach. The decomposition of the lattice into 
sublattices, which is the basis for the decomposition method, depends on the range of the interactions, so that 
this numerical approach is far less general than the predictor-corrector method. Crystal field anisotropies leave 
the performance of the predictor-corrector method almost unaffected, whereas the decomposition method suffers 
from a drastic reduction in speed. 

2. Advantages of the decomposition method 

The greatest advantage of the decomposition method is its capability for handling large time steps and the exact 
conservation of spin length. In the absence of anisotropies it also conserves the energy exactly and it maintains 
reversibility. For anisotropic Hamiltonians energy conservation and reversibility can be obtained to a high 
accuracy using iterative schemes; exact magnetization conservation, however, is lost. The time steps typically 
used are unaccessible by the predictor-corrector method due to the lack of exact energy conservation within 
the algorithm. The resulting speedup of a spin dynamics simulation gives access to larger lattices, provided 
a second order decomposition is sufficiently accurate for the problem under investigation. In simple cases the 
fourth order decomposition yields very accurate results even for time steps an order of magnitude larger than 
typical time steps used for the predictor-corrector method. 

A general recommendation for one or the other method cannot be given. If, however, special interactions are the 
main objective of the investigation, one should resort to the predictor-corrector method because of its generality. If the 
interactions in the system are standard two-spin exchange interactions and the objective is to study large systems, one 
should consider the decomposition method because of its speed and its built-in energy and spin length conservation. 
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Note added in proof 

While this work was in press we have learned that a method equivalent to the one presented here had been developed 
independently by Frank, Huang, and Leimkuhler (J. Frank, W. Huang, and B. Lcimkuhlcr, Journal of Computational 
Physics 133, 160 (1997).) Frank, Huang, and Lcimkuhlcr denote their method as Staggered Red-Black Method and 
employ the same decomposition technique which is described in this work. 
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FIG. 1. Energy e(t) — E(t) / (JL 3 ) per spin for the predictor-corrector method (see Eqs.(2.1) and ( |2.2j )) for the time step 
St = 0.01/J and D = (solid line). Both decomposition schemes (see Eqs.(3.5) and (3.6)) conserve the energy exactly (dashed 
line). 



Magnetization m(t) — \lSA(t)\/L 3 per spin for the predictor-corrector method (see Eqs.(^.l|) and (2.2)) for the time 

0.04/ J (solid 
5D yields fluctuations in 



FIG. 2. 

step St — 0.01/J (dashed line) and the second order decomposition scheme (see Eq.(3.5)) for the ti me s tep St 
line) for D — 0. The predictor-corrector method conserves the magnetization exactly, whereas Eq 
m(t) on all time scales. 



FIG. 3. Magnetization m(t) per spin for the fourth o rde r decomposition method (see Eq.(3.6)) with St — 0.2/J (solid line) 
and for the second order decomposition method (see Eq.(3.5)) with St — 0.04/J (dashed line) for D = 0. The CPU time needed 
for both methods with these time steps is almost the same. Despite the large time step the fourth order decomposition method 
yields substantially smaller magnetization fluctuations than the second order decomposition. 



FIG. 4. Energy e{t) per spin for the predictor-corrector method (see Eqs.( |2.l[ ) and ([2.2])) for the time step St = 0.01/J and 
D — J (solid line). The decomposition scheme given by Eq.(3.5) for St — 0.04/J shows a decrease of e(t) (dashed line) and 
Eq.(3.6)) for St = 0.2/J also yields a decrease of e(t) (dotted line). For these parameters the decomposition methods show 
better energy conservation than the predictor-corrector method. 



FIG. 5. Magnetization m z (t) per spin for the predictor-corrector method (see Eqs.( |2.l[ ) and (2.2)) for the time step 
St — 0.01/J (dashed line) and the second order decomposition scheme (see Eq.(p.q)) for the time step St = 0.04/J (solid line) 
and D = J. The qualitative behavior of m z (t) closely resembles the behavior in the isotropic case shown in Fig.B. 



FIG. 6. Magnetization m z (t) per spin for the fourth order decomposition method (see Eq.(3.6)) with St = 0.2/J (solid line) 
and for the second order decomposition method (see Eq.(3.5)) with St — 0.04/J (dashed line) for D — J. The typical temporal 
structure of the fluctuations of m z (t) is very similar to the isotropic case shown in Fig.^. 



FIG. 7. Longitudinal dynamic structure factor 5;(q, uj) of an isotropic Heisenberg ferromagnet for T = 0.8T C and |q| = n/5 
in the (100) direction on a simple cubic lattice (L = 10, periodic boundary conditions) normalized to the static structure factor 
Sj(q) = J °° Si(q,uj)du) for the predictor-corrector method (diamonds) and the second order decomposition method (triangles) 
for time steps St = 0.01/J and St = 0.04/J, respectively (see main text). The statistical error (one standard deviation) is 
smaller than the symbol size. 
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FIG. 8. Transverse dynamic structure factor £ t (q, ui) of an isotropic Heisenberg ferromagnet for T = 0.8T C and |q| = 7r/5 
in the (100) direction on a simple cubic lattice (L = 10, periodic boundary conditions) normalized to the static structure factor 
St(q) = <St(q, Lo)dw for the predictor-corrector method (diamonds) and the second order decomposition method (triangles) 
for time steps St = 0.01/ J and St = 0.04/ J, respectively (see main text). The error bars represent one standard deviation. 
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